// --- Includes and pre-defined constants ---
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <iostream>
#include <fstream>
#include <ctime>
#include <strings.h>
#include <cstring>
#include <cfloat>
#include <iomanip>
#include <algorithm>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <sys/time.h>
#include <unistd.h>
#define ERROR -1
#define ROW 1001
#define COL 1002
#define ALL 1003

using namespace std;

// --- Helper methods for double arrays ---


/*  Allocates memory for a 2-d array
 */
void allocate(double ***array, int rows, int cols){
	*array = new double* [rows]; // pointers to the rows
	(*array)[0] = new double [rows*cols]; // actual data block
	for (int i=1; i<rows; i++)
		(*array)[i] = (*array)[i-1] + cols; // point the row-pointers to appropriate data block locations
}

/*  Allocates memory for a 3-d array1
 */
void allocate(double ****array, int dim1, int dim2, int dim3){
	*array = new double** [dim1]; // pointer to the first dimension
	for (int i=0; i<dim1; i++)
		allocate(&((*array)[i]), dim2, dim3);
}

/*  Allocates memory for a 4-d array
 */
void allocate(double *****array, int dim1, int dim2, int dim3, int dim4){
	*array = new double*** [dim1]; // pointer to the first dimension
	for (int i=0; i<dim1; i++)
		allocate(&((*array)[i]), dim2, dim3, dim4);
}

/*  Allocates memory for a 1-d array
 */
void allocate(double **array, int length){
	*array = new double[length]; // allocate data block
}

/*  Frees up memory for a 2-d array
 */
void deallocate(double ***array){
	delete [] (*array)[0]; // free up block of doubles
	delete [] (*array); // free up row pointers
}

/*  Frees up memory for a 3-d array
 */
void deallocate(double ****array, int dim1){
	for (int i=0; i<dim1; i++)
		deallocate(&((*array)[i]));
	delete [] (*array);
}

/*  Frees up memory for a 4-d array
 */
void deallocate(double *****array, int dim1, int dim2){
	for (int i=0; i<dim1; i++)
		deallocate(&((*array)[i]), dim2);
	delete [] (*array);
}

/*  Frees up memory for a 1-d array
 */
void deallocate(double **array){
	delete [] (*array); // free up block of doubles
}

/*  Prints a 1-d array to standard output
 */
void print(double *array, int length){
	for (int i=0; i<length; i++)
		cout << array[i] << " " ;
	cout << endl;
}
void print(unsigned int *array, int length){
	for (int i=0; i<length; i++)
		cout << array[i] << " " ;
	cout << endl;
}

/*  Prints a 2-d array to standard output
 */
void print(double **array, int rows, int cols){
	for (int i=0; i<rows; i++){
		for (int j=0; j<cols; j++)
			cout << array[i][j] << " " ;
		cout << endl;
	}
}

/*  Prints transpose of a 2-d array to standard output
 */
void printTranspose(double **array, int rows, int cols){
	for (int i=0; i<cols; i++){
		for (int j=0; j<rows; j++)
			cout << array[j][i] << " " ;
		cout << endl;
	}
}

/*  Transposes a 2-d array
 */
void transpose(double **destination, double **source, int sourceRows, int sourceCols){
	for (int i=0; i<sourceRows; i++)
		for (int j=0; j<sourceCols; j++)
			destination[j][i] = source[i][j];
}

/*  Beta RNG for a 1-d array
 */
void betarnd(double *array, int length, double a, double b, gsl_rng* r){
	for (int i=0; i<length; i++){ 
		array[i] = gsl_ran_beta(r, a, b); 
	}
}

/*  Beta RNG for a 2-d array
 */
void betarnd(double **array, int rows, int cols, double a, double b, gsl_rng* r){
	for (int i=0; i<rows; i++){ 
		for (int j=0; j<cols; j++){
			array[i][j] = gsl_ran_beta(r, a, b);
		}
	}
}

/*  Gamma RNG for a 1-d array
 */
void gammarnd(double *array, int length, double a, double b, gsl_rng* r){
	for (int i=0; i<length; i++){ 
		array[i] = gsl_ran_gamma(r, a, b); 
	}
}

/*  Gamma RNG for a 2-d array
 */
void gammarnd(double **array, int rows, int cols, double a, double b, gsl_rng* r){
	for (int i=0; i<rows; i++){ 
		for (int j=0; j<cols; j++){
			array[i][j] = gsl_ran_gamma(r, a, b);
		}
	}
}

/*  Sum up a 1-d array
 */
double sum(double *array, int length){
	double answer=0;
	for (int i=0; i<length; i++)
		answer += array[i];
	return answer;
}

/* Sum up a 2-d array
 */
double sum(double *result, double **array, int rows, int cols, int direction){
	if (direction == ROW) {
		for (int i=0; i<rows; i++){
			result[i] = 0;
			for (int j=0; j<cols; j++)
				result[i] += array[i][j];
		}
		return log(0); // this return should not be used
	} else if (direction == COL) {
		for (int j=0; j<cols; j++)
			result[j] = 0;
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				result[j] += array[i][j];
		return log(0); // this return should not be used
	} else if (direction == ALL) {
		double answer=0;
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				answer += array[i][j];
		return answer;
	}
	return ERROR;
}

/* Sum up a 2-d array
 */
double sum(double *result, int **array, int rows, int cols, int direction){
	if (direction == ROW) {
		for (int i=0; i<rows; i++){
			result[i] = 0;
			for (int j=0; j<cols; j++)
				result[i] += array[i][j];
		}
		return log(0); // this return should not be used
	} else if (direction == COL) {
		for (int j=0; j<cols; j++)
			result[j] = 0;
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				result[j] += array[i][j];
		return log(0); // this return should not be used
	} else if (direction == ALL) {
		double answer=0;
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				answer += array[i][j];
		return answer;
	}
	return ERROR;
}

/*  Scale up a 1-d array by a constant
 */
void scale(double z, double *array, int length){
	for (int i=0; i<length; i++)
		array[i] *= z;
}

/*  Scale up a 2-d array by a constant
 */
void scale(double z, double **array, int rows, int cols){
	for (int i=0; i<rows; i++)
		for (int j=0; j<cols; j++)
			array[i][j] *= z;
}

/*  Scale up a 2-d array by a 1-d array
 */
void scale(double *z, double **array, int rows, int cols, int direction){
	if (direction == ROW) {
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				array[i][j] *= z[i];
	} else if (direction == COL) {
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				array[i][j] *= z[j];
	}
}

/*  Scale down a 2-d array by a 1-d array
 */
void scaleDown(double *z, double **array, int rows, int cols, int direction){
	if (direction == ROW) {
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				array[i][j] /= z[i];
	} else if (direction == COL) {
		for (int i=0; i<rows; i++)
			for (int j=0; j<cols; j++)
				array[i][j] /= z[j];
	}
}


/*  Add to a 1-d array
 */
void add(double z, double *array, int length){
	for (int i=0; i<length; i++)
		array[i] += z;
}

/*  Add to a 2-d array
 */
void add(double z, double **array, int rows, int cols){
	for (int i=0; i<rows; i++)
		for (int j=0; j<cols; j++)
			array[i][j] += z;
}

/*  Element-by-element Multiplication of two 1-d arrays
 */
void mult(double *result, double *array1, double *array2, int length){
	for (int i=0; i<length; i++)
		result[i] = array1[i]*array2[i];
}

/*  Element-by-element Multiplication of two 2-d arrays
 */
void mult(double **result, double **array1, double **array2, int rows, int cols){
	for (int i=0; i<rows; i++)
		for (int j=0; j<cols; j++)
			result[i][j] = array1[i][j] * array2[i][j];
}

/*  Dot product of two 1-d arrays
 */
double dotProduct(double* array1, double* array2, int length){
	double answer = 0;
	for (int i=0; i<length; i++)
		answer += array1[i]*array2[i];
	return answer;
}

/*  Dot product of two 1-d arrays
 */
double dotProduct(double* array1, int* array2, int length){
	double answer = 0;
	for (int i=0; i<length; i++)
		answer += array1[i]*array2[i];
	return answer;
}

/*  Copy for 3-d arrays
 */
void copy(double ***destination, double ***source, int dim1, int dim2, int dim3){
	for (int i=0; i<dim1; i++)
		for (int j=0; j<dim2; j++)
			for (int k=0; k<dim3; k++)
				destination[i][j][k] = source[i][j][k];
}

/*  Copy for 2-d arrays
 */
void copy(double **destination, double **source, int rows, int cols){
	for (int i=0; i<rows; i++)
		for (int j=0; j<cols; j++)
			destination[i][j] = source[i][j];
}

/*  Copy for 1-d arrays
 */
void copy(double *destination, double *source, int length){
	for (int i=0; i<length; i++)
		destination[i] = source[i];
}

/*  Set to a value for 1-d arrays
 */
void setValue(double z, double *array, int length){
	for (int i=0; i<length; i++)
		array[i] = z;
}

/*  Set to a value for 2-d arrays
 */
void setValue(double z, double **array, int rows, int cols){
	for (int i=0; i<rows; i++)
		for (int j=0; j<cols; j++)
			array[i][j] = z;
}

/*  Horizontal concatenation for 2-d arrays
 */
void horzcat(double **result, int rows, double **array1, int cols1, double **array2, int cols2){
	for (int i=0; i<rows; i++){
		for (int j=0; j<cols1; j++)
			result[i][j] = array1[i][j];
		for (int j=0; j<cols2; j++)
			result[i][cols1+j] = array2[i][j];
	}
}

/*  Vertical concatenation for 2-d arrays
 */
void vertcat(double **result, int cols, double **array1, int rows1, double **array2, int rows2){
	for (int i=0; i<rows1; i++)
		copy(result[i], array1[i], cols);
	for (int i=0; i<rows2; i++)
		copy(result[rows1+i], array2[i], cols);
}

/*  Computes time difference between two timeval structures
 */
double timeDiff(timeval endTime, timeval startTime){
	double answer = (double)(endTime.tv_sec) + endTime.tv_usec/1000000.0 - (double)(startTime.tv_sec) - startTime.tv_usec/1000000.0 ; 
	return answer;
}

/*  Add two numbers in log-space
 */
double logAdd(double a, double b){
	double answer = 0;
	if (a > b)
		answer = a + log(exp(b-a)+1);
	else if (b > a)
		answer = b + log(exp(a-b)+1);
	else
		answer = log(2)+a;
	return answer;
}

/*  My implementation of isinf()
 *  Returns -1 for -Inf
 *  Returns +1 for +Inf
 *  Returns 0 otherwise
 */
int myIsInf(double a){
	if (fpclassify(a) == FP_INFINITE){
		if (a<0)
			return -1; // -Inf
		else
			return +1; // +Inf
	} else {
		return 0; // finite double
	}
}  

/*  Minimum element in an array while ignoring -Inf/+Inf/Nan
 */
double* minIgnoreNegInf(double *array, int length){
	double min=DBL_MAX; double* answer=NULL; // start with min=max_double and pointer to null
	for (int i=0; i<length; i++){
		// if ((isinf(array[i]) != -1) && (array[i]<min)){
		if ((isfinite(array[i]) != 0) && (array[i]<min)){
			min = array[i];
			answer = &(array[i]);
		}
	}
	if (answer == NULL){
		cerr << "minIgnoreNegInf could not find min." << endl;
		print(array, length); cout.flush();
	}
	return answer;
}

/*  Add an array of numbers in log-space
 */
double logAdd(double *array, int length){
	double answer = 0;
	double * _array = new double[length]; copy(_array, array, length); // make a copy of the array
	double l_realmax = log(DBL_MAX) - log(length);
	double* minp = minIgnoreNegInf(_array, length); // find min ignoring -Inf
	double* maxp = max_element(_array, _array+length); // find max
	while (*maxp - *minp > l_realmax){ // check for very big difference between min and max
		*minp = log(0); // set current min to -Inf
		minp = minIgnoreNegInf(_array, length); // find the next min element
	}
	for (int i=0; i<length; i++){
		answer += exp(_array[i] - *minp); // add after removing common multiplier
	}
	answer = *minp + log(answer); // take in to account common multiplier
	delete [] _array;
	return answer;
}


/*  Computes mean over a 2-d array
 */
void mean(double* result, int start, int end, double** array, int rows, int cols){
	setValue(0, result, cols);
	for (int i=start; i<=end; i++)
		for (int j=0; j<cols; j++)
			result[j] += array[i][j];
	scale(1.0/(end-start+1), result, cols);
	
}

/*  Computes mean over a 3-d array
 */
void mean(double** result, int start, int end, double*** array, int dim1, int dim2, int dim3){
	setValue(0, result, dim2, dim3);
	for (int i=start; i<=end; i++)
		for (int j=0; j<dim2; j++)
			for (int l=0; l<dim3; l++)
				result[j][l] += array[i][j][l];
	scale(1.0/(end-start+1), result, dim2, dim3);
}

/* reads integer data from file */
int readData(int n, int d, char *inFileName, int** data){
	
	// open the data file for reading
	FILE* inFile = fopen(inFileName, "r"); // open the data file for reading
	if (inFile == NULL){
		cerr << "Could not open the data file " << inFileName << " for reading." << endl;
		return -1;
	} // error handling for opening the file
	
	// read from the data file
	int count = 0; // actual number of records read
	while (!feof(inFile) && count < n) { // read one line at a time
		char line[25000];
		fscanf(inFile, "%s", line); // read one line
		//cout << "read: " << line << endl;
		char* toks = strtok(line,",");
		// cout << "first token: " << toks << endl;
		for (int j=0; j<d; j++){
			toks = strtok(NULL,",");
			data[count][j] = atoi(toks);
		}
		//print(data[count],d);
		count++;
	}
	
	// close the data file
	fclose(inFile);	
	cerr << "Number of data points ACTUALLY read: " << count << endl;
	return 0;
	
}

int readMu(int k, int d, char *inFileName, double** mu){
	
	// open the data file for reading
	FILE* inFile = fopen(inFileName, "r"); // open the data file for reading
	if (inFile == NULL){
		cerr << "Could not open the mu file " << inFileName << " for reading." << endl;
		return -1;
	} // error handling for opening the file
	
	// read from the data file
	int count = 0; // actual number of records read
	while (!feof(inFile) && count < k) { // read one line at a time
		char line[25000];
		fscanf(inFile, "%s", line); // read one line
		//cout << "read: " << line << endl;
		char* toks = strtok(line,",");
		mu[count][0] = atof(toks); 
		for (int j=1; j<d; j++){
			toks = strtok(NULL,",");
			mu[j][count] = atof(toks);
		}
		//print(data[count],d);
		count++;
	}
	
	// close the data file
	fclose(inFile);	
	cerr << "Number of mu points ACTUALLY read: " << count << endl;
	return 0;
	
}

/*  Prints a 2-d array to file
 */
void print(double **array, int rows, int cols, char* outFileName){
	ofstream outfile;
	outfile.open(outFileName);
	for (int i=0; i<rows; i++){
		for (int j=0; j<cols; j++)
			outfile << array[i][j] << " " ;
		outfile << endl;
	}
	outfile.close();
}

/*  Prints a 1-d array to file
 */
void print(double *array, int length, char* outFileName){
	ofstream outfile;
	outfile.open(outFileName);
	for (int i=0; i<length; i++)
		outfile << array[i] << " " ;
	outfile << endl;
	outfile.close();
}

/* reads integer data with time stamps as first token */
int readDataWithTime(int n, int d, char *inFileName, int** data, int* time){
	
	// open the data file for reading
	FILE* inFile = fopen(inFileName, "r"); // open the data file for reading
	if (inFile == NULL){
		cerr << "Could not open the data file " << inFileName << " for reading." << endl;
		return -1;
	} // error handling for opening the file
	
	// read from the data file
	int count = 0; // actual number of records read
	while (!feof(inFile) && count < n) { // read one line at a time
		char line[25000];
		fscanf(inFile, "%s", line); // read one line
		//cout << "read: " << line << endl;
		char* toks = strtok(line,",");
		time[count] = atoi(toks);
		// cout << "first token: " << toks << endl;
		for (int j=0; j<d; j++){
			toks = strtok(NULL,",");
			data[count][j] = atoi(toks);
		}
		//print(data[count],d);
		count++;
	}
	
	// close the data file
	fclose(inFile);	
	cerr << "Number of data points ACTUALLY read: " << count << endl;
	return 0;
	
}

/* reads double data with time stamps as first token */
int readDataWithTime(int n, int d, char *inFileName, double** data, int* time){
	
	// open the data file for reading
	FILE* inFile = fopen(inFileName, "r"); // open the data file for reading
	if (inFile == NULL){
		cerr << "Could not open the data file " << inFileName << " for reading." << endl;
		return -1;
	} // error handling for opening the file
	
	// read from the data file
	int count = 0; // actual number of records read
	while (!feof(inFile) && count < n) { // read one line at a time
		char line[25000];
		fscanf(inFile, "%s", line); // read one line
		//cout << "read: " << line << endl;
		char* toks = strtok(line,",");
		time[count] = atoi(toks);
		// cout << "first token: " << toks << endl;
		for (int j=0; j<d; j++){
			toks = strtok(NULL,",");
			data[count][j] = atof(toks);
		}
		//print(data[count],d);
		count++;
	}
	
	// close the data file
	fclose(inFile);	
	cerr << "Number of data points ACTUALLY read: " << count << endl;
	return 0;
	
}
